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ABSTRACT 

We discuss a new criterion to estimate the mass in the outer, non-equilibrium region of galaxy 
clusters, where the galaxy dynamics is dominated by an overall infall motion towards the 
cluster centre. In the framework of the spherical infall model the local mean velocity of the 
infalling galaxies at every radius provides information about the integrated matter overdensity 
5. Thus, a well-defined value of the overdensity 5t is expected at the turnaround radius r^, 
i.e. the radius where the Hubble flow balances the infall motion. Within this scenario, we 
analysed the kinematical properties of a large catalogue of simulated clusters, using both 
dark matter particles and member galaxies as tracer of the infall motion. We also compared 
the simulation with analytical calculation performed in the spherical infall approximation, to 
analyze the dependence of the results on cosmology in spatially flat universe. If we normalize 
cluster mass profiles by means of the turnaround mass Alt (i-S- the mass within r^), they are 
consistent with an exponential profile in the whole non-equilibrium region (0.5 < r/rt < 2). 
Turnaround radii are proportional to virialization radii (rt ~ 3.5ru), while turnaround masses 
are proportional to virialization masses, i.e. Mt ~ 1.7M„, where My is the mass within r^. 
Actually, the mass evaluated within the turnaround radius is a more exhaustive evaluation of 
the total mass of the cluster These results can be applied to the analysis of observed clusters. 

Key words: cosmology: large scale structure - galaxies: clusters: general - galaxies: kine- 
matics and dynamics 



1 INTRODUCTION 

The gravitational collapse of galaxies towards the centre of clusters 
is usually described within the framework of the spherical infall 
model, as the motion of a set of concentric, spherically symmetri- 
cal mass shells jsee e.g. GuTm & Gott 1972, Silk 1 97j)|Schechter| 
|1980| l. Actually the spherical infall model is widely accepted in 
literature, since it describes fairly well the dynamics of the non- 
equilibrium region of galaxy clusters, defined as the region where 
the effects of virialization and the crossing of the above-mentioned 
shells are negligible and some overall infall motion of member 
galaxies is recognizable. Under the spherical symmetry assump- 
tion, the infall motion produces a pattern of caustic surfaces in the 
galaxy redshift-space distribution (which is obtained representing 
the line-of-sight velocities of galaxies versus their projected posi- 
tion on the sky plane). These caustics envelop all galaxies whose 
infall motion overwhelms the Hubble flow ( |Kaiser|1987| l. Caustics 
with a characteristical "trumpet" shape were actually observed in 
the redshift-space distribution of clusters (see e.g. Ostriker et al.| 
[1998). Diaferio & Geller piaferio & Geller„1997) and Diaferio 
I Diaferio|1999^ showed that the caustic amplitude provides a direct 
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measure of the escape velocity of galaxies, and therefore allows to 
estimate the mass profile of the cluster in the innermost part of the 
non-equilibrium region, up to the turnaround radius rt (i.e. the ra- 
dial distance where the velocity of the infall motion is equal to the 
Hubble flow velocity.). The caustic technique was applied to the ob- 
servation of many local clusters ( jsee e.g. Geller, Diaferio, & Kurtz| 
1999'; "Rmes & al. 2000, 2001a, 20031. These mass estimates are 
consistent with those based on virial theorem ( Girardi et al. 11998} 
Biviano & Girardi 200^ and weak tensing observations ( [Diaferio^ 
|Geller, & Rines|2006j and references therein). In fact, up to now 
the sampled volumes were always restricted within the turnaround 
radius, and this is due to the definition of caustics surfaces ( |Reg6s| 
|& Geller] 1989l l. 

In this paper, we discuss an approach to the issue of mass es- 
timation, which can be applied to larger sampled volumes, well be- 
yond the turnaround radius. We use the radial velocity of galaxies 
as the key quantity, instead of the escape velocity, as in the caustic 
technique. According to Silk (|Silk 1974), Peebles ( |PeebIes|19"76| 
[T980) , and Gunn ( |Gunn|1978^ , within the spherical symmetry hy- 
pothesis, the velocity of the matter infall motion at a certain dis- 
tance from the centre depends on the encompassed mass. Our pur- 
pose is to use this dependence to constrain the value of the over- 
density at the turnaround radius. In fact, we see that the turnaround 
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radius is far outside tiie virialization core of clusters, and is thiere- 
fore a suitable normalization scale for the cluster mass profile in 
the non-equilibrium region ( |Vedel & Hartwick|1998 i To test our 
assumptions and to verify the results, we will analys a large galaxy 
population extracted from a simulated cluster catalogue (Borgani 
|et al.|2004||Biviano et al.|2006| l. We will study all clusters both as 
a whole and one by one. We will prove that the actual turnaround 
overdensity of clusters is in good agreement with the predictions 
of the spherical infall model. Moreover, we will show that the nor- 
malized mass profiles are generally consistent with a power-law 
profile, which extends the standard Navarro-Frenk-White profile 
( |Navarro, Frenk, & White 1995] [|[l996] [1997] hereafter NEW) to 
the non-equilibrium region. 

In Section [2] we present the details of our model, concerning 
the theoretical framework l |2.1| > and the simulated data sample \2.2) . 
In Section[3]we discuss the results of our analysis, focusing on the 
mass estimation at the turnaround radius \2>A) and in the whole 
non-equilibrium region, up to 8 virialization radii l |3.2^ . Finally, in 
Section|4|we draw the conclusions of our work. 



2 THE MODEL 

2.1 Theoretical framework 

Consider a galaxy located at a distance r from the centre of a clus- 
ter We call 'infall velocity' Vr the peculiar velocity of the galaxy 
along the radial direction (i.e. towards the cluster centre), assuming 
that it is positive when directed inwards. The matter overdensity 5 
is defined with respect to the background density pjg = ^opcr as 
follows: 
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where fio is the cosmological matter density parameter and per is 
the critical density; all the quantities are considered at the present 
day. According to the hypotheses of the spherical infall model, the 
ratio between the infall velocity and the Hubble flow velocity Har 
(where Hq is the Hubble parameter) can be written unambiguously 
as a function F of both fio and 5 ( |Silk| 1 974[ [Peebles| 1 976[ [Gurinj 
[T978j|Peebles|1980| l: 
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Several definition of F were proposed in literature. Regos and 
Geller jRegos & Geller|1989) demonstrated that quite for all pur- 
poses F may be factored into a polynomial P of the mere overden- 
sity 5: F{Q.o, 5) ~ Q.l-^P{S). Lightman & Schechter i Lightman & 



|Schechter|1990^ described a simple approach to compute the high- 
order polynomial terms. However, this formulation is devised to fit 
the spherical infall in the very-low-overdensity region (5 < 2). A 
better agreement in the whole non-equilibrium region {5 < 30) 
is obtained with non-polynomial approximations ( YahilT985' 'Vil-| 
lumsen & Davis|1986| see Section [TT) . Lahav et al. (Lahav et al. 



1991| l took into account the possibility of a non-zero cosmological 



constant parameter Aq at the present day, and obtained a corrective 
term accounting for a 3-per-cent discrepancy with the previous re- 
sults. Due to the small size of this correction, we will not consider 
here the effect of Aq, and will assume hereafter that F is approx- 
imately factorable into a cosmological term Vl'l'^ and a generical 
function / of 5: 
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Equation ^ and equation ^ were commonly used to evalu- 
ate Qo from observations of local clusters i Regos & Geller 1989] 
|Lynde n-BeU, Lahav, & Burstein||1989l|Lahav et al.||1991|l. Vice 



versa, since we are handling a simulation and therefore we do 
know the cosmology, we can reverse this approach and use the 
equations to compute 5 as a function of Vr/ Hor. In principle, in 
a purely spherically-symmetric scenario, this would constrain the 
whole overdensity profile along r, because there would be a one- 
to-one dependence between r and Vr/Hor. But this assumption 
is actually too restrictive to describe the overall infall motion of 
galaxies, even in the non-equilibrium region. In fact, the presence 
of small-scale substructure is proved to locally affect the galaxy 
velocity, introducing a sort of random "kinematical" noise which 
blurs the infall velocity profile i jDiaferio & Geller||1997[ l. Never- 
theless, we will prove that this profile is regular enough to make 
possible the estimate of the turnaround radius rt, which is defined 
by the condition Vr/Hort ~ 1. If so, equation ^ can be used to 
implicitly define the turnaround overdensity St = S{rt) as a func- 
tion of Qq: 
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Since we are using the cosmological N-body simulation de- 
scribed in Section [2^ we have analyzed the dependence on flo by 
means of analytical calculations, which assume spherically sym- 
metric infall, and a spatially flat Universe. This model provides an 
analytical definition of Vt and 5t for a cluster jEke et al.|1996|[La^ 
|havetal.|199"T] see App.|A|: 
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In these equations, r„ is the virialization radius of the cluster, while 
the constants and Kt parametrize the amplitude of two density 
perturbations which are respectively collapsing and turning around 
at the present day. The angle Ov and the parameter are both re- 
lated to K„, while the angle 6t is related to nt (see App.|A]for fur- 
ther details). In principle, all parameters in equation ^ and equa- 
tion |6]( depend on the adopted cosmology, and in particular on the 
value of Qo - In the range 0.2 5C fio ^ 0.4, the dependence on fio 
is approximately factorable as follows (with a one-per-cent accu- 
racy): 
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where f t and St are independent from cosmology. According to the 
latest WMAP observations fSpe rgel et al.|2003^ , Qo = 0.27±0.04. 
Substituting this value into equation Q and equation (|6|, we get 
rt = (3.02 ± 0.02)r„ and St = 12.2l^ j, corresponding to 
rt = (0.72 ± 0.03)ft and St = (2.7 ± 0.3)5t from equation ^ 
and equation ijsjl, where ft = 4.2 and St = 4.6. As one can see, the 
uncertainty on rt and St due to dependence on cosmology is very 
small and can be neglected for our purposes (as it will be shown in 
Section ing . For these reason, we will adopt hereinafter the con- 
cordance value Qo — 0.3 if not otherwise specified. 

The turnaround radius rt will be adopted as a normalization 
scale useful to describe the outskirts of clusters; in this way, it 
replaces the virialization radius, which is usually adopted in the 
cluster core. Within this framework, we will demonstrate that in 
the non-equilibrium region the overdensity profile (5ne('") and the 
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Figure 1. Characteristics of the sample. First plot: frequency distribution of 
the virialization radii r„ ; second plot: frequency distribution of the virial- 
ization masses My. n^-i^ci is the number of cluster per frequency bin. 



mass profile Mne('') of a given cluster are generally consistent 
with a single profile, if they are normalized to the turnaround scale. 
Therefore, we can write: 



5NE(r) = (l + St) 



rt 



ffNE(r) - 1, (9) 



MNE(r) = MtQNEir) 



4 3 

-■KTt r2oPcr(l + 5t)ffNE(r-), (10) 

where (;ne is a function to be defined. Equation ([9]l and equation 
formally correspond to the equations of the Navarro-Frenk- 
White profile l |Navarro, Frenk, & WhitefT995l||1996||T997| here- 
after NEW). The difference between our profile and the NEW one 
lies in the choice of the normalization scale (the turnaround radius 
rt and the turnaround overdensity 5t substitute the virialization ra- 
dius r„ and the virialization overdensity 5^), and lies also in the 
definition of ^ne, which we will see to be different from the corre- 
sponding NEW function (7nfw('') = ln(l + Ci,r/r„) — Ci,r/(ci,r + 
r„), where c„ is the cluster concentration parameter i see e.g. Bul- 
|lock et al.|200T||Lokas & Mamon|200"T] i. 

2.2 The simulated catalogue 

The model was tested on a catalogue of 1 14 simulated clusters, with 
an overall population of 9631 galaxies. The clusters and the galax- 
ies were extracted by Biviano et al. ^Biviano et al.||200 6l from a 
large cosmological hydrodynamical simulation performed by Bor- 
gani et al. ( [Borgani et al.|2004) . We refer to these papers for a de- 
tailed description of the data sample. We just remark what follows: 

(i) The simulation was run with the tree-SPH GADGET-2 code 
( |Springel, Yoshida, & White ' 20011 fSpringel & Hemquist||2002l l, 
adopting a A-CDM cosmology (f^o = 0.3, Q.a = 1 — f2o, 
^bar = 0.019ft"^, h = 0.7 and erg = 0.8). It traced the evo- 
lution of 480"^ dark matter (DM) particles and 480'' gas particles 
(partly converted into stellar particles during the run) within a box 
of volume (192/1"^)^ Mpc^ 



(ii) The clusters were identified at 2; = with a standard 
Eriends-of-Eriends (EoE) algorithm, taking into account the DM 
particles of the simulation. After the identification, a spherical over- 
density algorithm was applied to determine the size of the virializa- 
tion core of each cluster. The virialization overdensity was defined 
as follows, in agreement with the adopted cosmology jBryan &] 
|Norman|1998l l: 

(1 + 5y)no = IStt^ + 82(^0 - 1) - 39(!^o - 1)^ ^ 101. (11) 

(iii) The galaxies were identified with the publicly available al- 
gorithm SKID <S tadel|2001| l; in this case, only the stellar component 
was taken into account. 

According to the definition of (5„ , we will define the virializa- 
tion radius and the virialization mass of each cluster as r„ = rioi 
and Mv = i\/ioi, respectively. The extracted clusters are very dif- 
ferent in size, with ranging from 0.88/i^^ Mpc to 2.23ft^^ Mpc 
and M„ ranging from 7.95 x 10^^/i~^Mq to 1.30 x 10^^/i^^Mq. 
The frequency distribution of r„ and M„ among the sample is 
shown in Eig. [T| Since many authors prefer to use r2oo and M200 
instead of and i\f„, we provide the average ratios r„/r2oo and 
Af„/M2oo computed on the entire cluster catalogue: 



1.36 ±0.04, 



(12) 



1.26 ±0.11. 



(13) 



'"200 

A/. 
A/200 

The number of member galaxies is very different in different 
cluster, ranging from 17 to 403. To make different objects compa- 
rable, we sliced all clusters into a set of concentric shells, using the 
virialization radius as the scale reference. The shells were defined 
so as to cover the whole extent from the virialization core to the far 
outskirts of clusters (i.e. from 0.1r„ to 8r„). We adopted a logarith- 
mical spacing in order to fit the decreasing galaxy number density 
along the radial coordinate. We defined the outer radius rj of each 
shell j as follows: 

-,{j/50)-l 



J = 1, 



,91. 



(14) 



(The choice of 91 shells is technical, induced by the quality of our 
data.) The same spacing was used to reconstruct the matter distri- 
bution along the radial coordinate in the clusters. 

The integrated overdensity and mass profile were computed 
taking into account all the particles in the simulation (i.e., DM, gas, 
and stellar particles), while the infall velocity was extracted from 
the member galaxies alone. We adopted this approach to better in- 
vestigate the possibility of reconstructing the cluster mass distribu- 
tion using only the dynamical properties of the member galaxies, 
which at least in principle can be directly inferred from the obser- 
vations. Within this approach, 5j and Mj are defined respectively 
as the overdensity and the mass of all the particles enclosed within 
the sphere of radius r-j, while Vr;j is the mean infall velocity of 
all the galaxies within the shell j. We used the DM particles as a 
tracer of the cluster dynamics only when computing the 3-d veloc- 
ity dispersion within the virialization core of clusters (Section [3T| , 
because in this case the DM component yields a stabler result due 
to its larger statistical significance. 



3 THE RESULTS 

3.1 The turnaround radius and the overdensity estimation 

Eig. [2] and Eig. [3] represent the sample distribution of the normal- 
ized infall velocity Vr/ Hor of member galaxies as a function of the 
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Figure 2. Overall infall velocity profile of member galaxies as a function 
of the normalized radial distance. The distribution of galaxies (points) has 
been smoothed with a running median (thick solid line) and interpolated 
with a power law (dashed line; see text). The horizontal bar indicates the 
turnaround condition tir/^^o'' = 1- 
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Figure 3. Overall infall velocity profile of member galaxies as a function 
of the overdensity. The distribution of galaxies (points) has been sinoothed 
with a running median (thick solid line) and then interpolated with fu^ 
(dotted line), /y (dashed line), and /m (narrow solid line). The horizontal 
bar indicates the turnaround condition Vr/Hgr = 1. 



normalized radial distance r / r „ and the overdensity 5, respectively. 
We superimposed all clusters into a single synthetic object, in or- 
der to increase the statistical significance, as suggested by jVedel^ 
|Hartwick| ( | 1 998 1 . Each point corresponds to a single galaxy, while 
the thick solid line is the running median (RM) of the distribution. 
Although the large variance, we can well recognize a common pro- 



file in the intervals r > 2.2r-„ and S < 35. This trend indicates the 
existence of an overall and well-defined galaxy infall motion in the 
non-equilibrium region of clusters. 

In both Fig.[2]and Fig.|3] we use the RM as a reference profile 
to describe the overall dynamics in the non-equilibrium region, in 
order to find the values of the turnaround radius and the turnaround 
overdensity. The radius rt is computed by interpolating the RM 
profile of Fig. |2] with a power law (dashed line). A linear-fitting 
algorithm applied to the bilogarithmic distribution gives: 



logi 



Vr 

Hor 
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where a^^ = 0.96 ± 0.03 and /3„^ = 1.72 ± 0.02 (1-cr uncer- 
tainties). Equation l |15[ ) is in good agreement with the RM profile 
for r > 2.2r-„, and can be used to compute the median turnaround 
radius rt,RM, defined by the condition / Hort\j^yi = 1 (corre- 
sponding to the horizontal bar in Fig.|2j. We obtain: 



rt.RM = (3.61±0.02)r-„. 



(16) 



In this equation, the 1-cr uncertainty is due to the interpolation algo- 
rithm and does not take into account the variance among the clus- 
ters in the sample, which will be considered later on. 

To determine St, we compare the RM profile in Fig. [3] with 
the plots of three different expressions of the function / defined in 
equation |3](, namely the linear approximation fun ( |Peebles|1976[ 
[Gunn 1978, Peebles 1980 dotted line), the Yahil approximation 
/y lYahil 1985 dashed line), and the Meiksin approximation /m 
( fViilum sen & Davis|1986| narrow solid line): 
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According to equation (j4|, these expressions provide as many im- 
plicit definition of the turnaround overdensity, corresponding to the 
points of intersection, in Fig. [3] of the three curves with the hor- 
izontal bar Vr/Hor = 1. As one can see from Fig. [5] the lin- 
ear approximation poorly describes the infall motion in the non- 
equilibrium region, since it departs from the data distribution even 
in the low overdensity region. Conversely, the non-linear functions 
/y and /m are close to RM up to the turnaround region. Equation 
l |17| ( and equation ^19\ can be inverted by trivial algebraic compu- 
tation, while equation l |18[ ) requires an ad hoc treatment (see App. 
[B). The results obtained in the three cases are, respectively. 
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It is worth noticing that the Yahil approximation agrees with the 
prediction of the spherical infall model (see Section |2.1| l. Con- 
versely, the Meiksin approximation shows the best agreement with 
the simulated data, being very close to RM profile in the inter- 
val (5 < 40. This is an evidence of the disagreement between the 
purely-spherical description and the actual galaxy dynamics in the 
non-equilibrium region of galaxy cluster. 

We also estimated the value of rt and St for single clusters 
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Figure 4. Turnaround radius estimation for all clusters in the data sample. 
The individual value of r^.j extracted from the DM distribution of galaxies 
(crosses and en'or bars) is compared with rt^RM (dashed line). The eiTor 
bars (1 cr) correspond to the uncertainty associated to the fit of the infall 
velocity profile, i is the index number of each cluster. 
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Figure 5. Turnaround overdensity estimation for all clusters in the data sam- 
ple. The value of Sfj extracted from DM distribution (crosses and error 
bars) is compared with St u^ (dotted line), 5t,Y (dashed line), and 5t.M 
(solid line). The error bars (1 a) correspond to the uncertainty associated to 
the fit of the infall velocity profile, i is the index number of each cluster. 



separately. In this case, we traced the infall velocity profile using 
the DM particles instead of the galaxies, since the galaxy distribu- 
tion is typically noisier (because in a single cluster the number of 
galaxies is much smaller than the number of DM particles traced 
by the simulation). Fig.|4]and Fig. |5] represents the individual val- 
ues rt:i and Sf.i, respectively, as a function of the index number 



i of our cluster catalogue (crosses and error bars). The error bars 
(1 cr) correspond to the uncertainty associated to the fit of the in- 
fall velocity profile. Most of values lie in a quite narrow band. In 
particular, the turnaround radius is (in almost all cases) not only 
larger but considerably larger than the virialization radius, so it lies 
in the infall region, and therefore it can be adopted as a suitable 
normalization scale for the non-equilibrium region. Averaging the 
individual values over the whole sample we obtain: 



logio ( ^ 1 =0.54 ±0.05, 



logio (1 + St) = 1.2 ±0.1. 



(23) 



(24) 



corresponding to ft — (3.5 ± 0.4)ri, and 5t = I5I3. uncer- 
tainty on these estimates is much larger than the uncertainty related 
to the value of (see Section [Z7) ; for this reason the latter has 
been neglected, and the concordance value Qq = 0.3 has been 
used throughout. There is a highly significant correlation between 
logio and logjg(l ± St), as indicated by the Pearson's correla- 
tion coefficient rp — —0.61 (significance ^ 99%). This result is 
not surprising, since a mutual dependence between v,- and r, and 
between v,. and 5, yields naturally a mutual dependence between 
rt and St . Taking into account this correlation, we can estimate the 
turnaround mass Mt = M4rt/rtf{l + <5t)/(l + Sv) with the 
correct l-cr uncertainty, as follows: 



logio ( ^ ) =0.24 ±0.01, 



(25) 



corresponding to Mt = (1.74 ± 0.04)M„. Equation 1 25 1 confirms 



the estimate of |Rines & Diaferiol ( |2006| >, which measured the aver- 
age ratio Mt /M200 by analysing a sample of observed clusters. 

The value of ft from equation l |23[ > is in agreement with the 
value of rt, RM from equation l |16[ ) (dashed line in Fig.|4|, and the 
1-(T uncertainty on the turnaround radius now takes into account the 
variance among the clusters. On the other hand, equation \24\ rules 
out the linear approximation result 5t,iin (dotted line in Fig.Bjl, and 
confirms the non-linear estimates 5t,Y and St,M (dashed and solid 
line in Fig.|5] respectively). The Meiksin approximation provides 
again the best estimate; therefore, it will be adopted hereafter as the 
best expression for /. 

The values in equation l[23j and l |24[ ( are not different within 
the uncertainties from the corresponding predictions of the spheri- 
cal infall model. However, the marginal evidence (within 1 a) of 
larger value of ft could be ascribed to a collapse which is not 
perfectly spherical. In fact, Hoffman (1986) found that a shear in 
the velocity field can induce higher infall velocities, which lead to 
larger turnaround radii. 

Biviano et al. ( [Biviano et al.|[2006| > show that the virializa- 
tion mass Mv depends on the 3-dimensional velocity dispersion of 
the DM component within the virialization radius, CTu.dm- Since 
the ratio Mt/M^ is quite constant in our catalogue, we expect to 
find a relation between the turnaround mass and cr^, dm- We there- 
fore compare the estimated values of cluster turnaround masses 
Mt-i with the values of CTt,,DM extracted from the simulated cat- 
alogue. The dependence of cluster turnaround masses on the veloc- 
ity dispersion is shown in Fig.|6j the solid line is obtained applying 
a linear-fitting algorithm to the bilogarithmic distribution, which 
gives 



logi 



Ml 



,(26) 



=10 \ ioi4/i-iM„ ) ' """i" V 103km s-i - 

where ctM^ = 0.6 ± 0.1 and Pm^ = 2.4 ± 1.3. We also computed 
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Figure 6. Dependence of the estimated turnaround mass on the 3-d DM 
velocity dispersion within the virialization core. The crosses represent the 
expected values (cj, Mt;i), while the error bars represent the l-cr uncer- 
tainty on the estimation of Mf-i, where i is the index number of our cluster 
catalogue. These values are compared with the best-fitting power law (solid 
line) and the best-fitting cubic relation (dashed line). The cubic relation 
which best fits the corresponding distribution of viiialization masses M„ 
(Biviano et al. 2006) is represented by the dotted line. 



the best-fitting cubic relation for the same distribution (dashed 
line): 



logi 



Mt 



QMt + 3 logi 



lO^km s~ 



(27) 



where aMt = 0.73 ± 0.05. Equation {11) is consistent with equa- 
tion p.6\ within the uncertainties. The cubic relation is favoured by 
Biviano et al. i jBiviano et al.|20()6) to describe the Mv-ctv depen- 
dence (since M„ ~ cr„r„ and r„ ~ a„), and therefore it is ex- 
pected to work also to describe the Mt-(Jv dependence. Our value 
of OLMt is consistent with the corresponding value by Biviano et al. 
( [Biviano et al.[20 06 1. We point out that in principle it is possible to 
use equation \26) or equation l[27j to obtain a mass estimate entirely 
based on the 3-dimensional velocity dispersion of the DM particles. 
One could also use galaxies instead of dark matter, provided that 
the galaxy velocity dispersion is an unbiased estimator of cr„, dm- 
This point is still debated in the literature, see, e.g.,|Biviano et al. 



(2006) and references therein. However, in their analysis, Biviano 
[et al. ( 2006) found that the bias is negligible when all galaxies (not 
only early-type galaxies) are considered. 



3.2 Overdensity and mass profile estimation 

Once rt and 5t are known, we use them to normalize the overden- 
sity profile and the mass profile of the clusters. The normalized 
profiles 5 and M are obtained in two different ways (i is the label 
of the cluster and j is the label of the cluster shell considered; see 
Section|3TT): 

(1) Using the individual values of turnaround radius and 
turnaround overdensity, computed for each cluster separately: 




Figure 7. Radial profile of the normalized overdensity (5'^' and of the nor- 
malized mass A/(^) for the whole sample. The profiles were reconstructed 
through the distribution of member galaxies (points) and subsequently in- 
terpolated with 5iME and A/ne, respectively (solid lines); see text and equa- 
tions |33} and |34j. In both plots, we highlighted a regular profile (naiTow 
solid line), and two irregular profiles (dashed line). 



rty. 



1 + Si, J 

1 + 5f, ' 



M. 



(1) 



Mf, 



(28) 



where Mt;i = M^{rt;i/r^f{l + 5t;i)/{l + 5^); 

(2) Using the mean values of turnaround radius and turnaround 
overdensity, obtained from the data of the whole cluster sample: 



=(2) 



1 + StA 



Mi. 



M, 



(29) 



where Mt,M = My{rt.RM/ry)^{l + (5t,M)/(l + Sy). 
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Table 1. Value of the parameters in equation j30j for different choices of 
Co- 

= 0.20 Ho = 0.27 Qo = 0.30 Ho = 0.4 

u 5.49 5.51 5.52 5.53 

7 0.49 0.51 0.51 0.52 

e 0.89 0.86 0.85 0.82 



Table 2. Agreement with the model for different types of mass profiles. We 
used the condition in equation \?>5\ with the normalization in equation ^28^ . 

e = 0.10 e = 0.15 e = 0.20 

Regular 48/72(67%) 70/72 (97%) 72/72(100%) 
IiTegular 19/42(45%) 29/42 (69%) 38/42(91%) 
All 67/114(59%) 99/114(87%) 110/114(97%) 



These two normalization criteria are useful to test the reliability of 
our model when applied to a sample of imulated clusters, as we will 
discuss later. 

To compare the extracted profiles with our model, we need an 
explicit expression of the function (;ne. A possible expression is 
obtained via the spherical infall model, considering the evolution 
of a spherical perturbation from a primordal time ti to the present 
time to. According to |Lahav et aL] ( |1991| l and |Lilje & Lahav| ( [T99T] l, 
the overdensity profile of a primordial perturbation can be written 
in the following way: 



5{r,U) = 



D(U 



2-K^aor D(to) jq 
V — 7^1^ — 76* 



X 



1-72 



+ 



fcji(fcr)P(fc)e-<^^'='''''' X 
37(1-72) 



dk. 



(30) 



In this equation, P{k) is the power spectrum of the perturbation 
(Bardeen et al. 1986), D{t) is the growing solution of density fluc- 
tuations fCarroU, P ress, & Tumer|1992^ , ji is the first-order spheri- 
cal Bessel function, i?/ is a filtering scale, (Tq is the rms fluctuation 
of the filtered density field, and the parameters R,, i^, 7, and 6 are 
related to the number density of peaks in the filtered density field 
( |Lilje & Lahav|[r991^ . We refer to above quoted papers for fur- 
ther details. All the parameters were tuned to match the conditions 
of the simulation we used. We chose in particular erg = 0.8 and 
Rf = Q.lh~^ Mpc, to avoid smoothing of fluctuations at the Mpc 
scale. The expectation values of i/, 7, and 6 were calculated in or- 
der to reproduce the number density of clusters in the simulation 
(which adopts Q.(, — 0.3); moreover, the values reported in Table[T| 
evidence the pure dependence on Slo of u, 7, and 6. 

Once the primordial overdensity profile is known, the spher- 
ical infall model provides a way to compute the corresponding 
present-day profile. The details are discussed in App. |A] We ob- 
tain, in a comoving framework: 



1 + S(r, to) 



0.0 



{l + 5{r,U)), 



(31) 



where the ratio r/a at any time is obtained by numerical integration 
of the Friedmann equation for the perturbation. The function (;ne 
is thus determined as the best fit function for the present-day mass 
profile, according to equation ^10\ . We performed the computation 
for different values of f2o in the range 0.2 ^ ilo ^ 0.4. Thus, the 
best fit function we adopted is: 



5ne('") = exp 



K 



1/4 



^-1 



(32) 



where K = 0.6 ± 0.1 (the uncertainty is due to the fitting 
algorithm). Once again the dependence on cosmology is very 
small. Assuming the concordance value Qo we obtain gNE{r) — 
exp[0.8{r/rt — 1)]. This relation is adopted hereinafter to com- 
pute the theoretical profiles 5ne(t') and A'/ne('') from equation |9| 
and equation ijTOjl. We can also compute the normalized theoretical 
profiles as follows: 



5NE(r) =f~^exp[0.8(f- 1)], 
MNE(f) = exp[0.8(f - 1)]. 



(33) 
(34) 



The two plots of Fig.|7]show the comparison between the nor- 
malized profiles extracted from the simulated cluster sample and 
the corresponding normalized profiles predicted by our spherical 



model. As one can see, the distribution of both V and M) ■' is in 
good agreement with ^neC*^) and AfNE('^) (solid lines). The simu- 
lated profiles show an intrinsic variance, due to the possible pres- 
ence of different cluster substructure (filaments, clumps, or even a 
bimodal core). We can roughly distinguish two basic patterns: 

(i) Regular profiles, i.e. smooth monotonic profiles, 

(ii) Irregular profiles, i.e. profiles showing one or more impor- 
tant changes of slope. 

A typical regular profile and two irregular profiles are shown in 
Fig. |7] with narrow solid lines and dashed lines, respectively. We 
qualitatively recognized in our sample 72 regular profiles and 42 
irregular profiles (about 63% and 37%, respectively). Usually, reg- 
ular profiles are expected to produce the best agreement with our 
integrated model, while profiles with changes of slope are expected 
to deviate from our prediction. 

We tested the behaviour of the profiles by analysing the pro- 
files of all the clusters of our catalogue when normalized by the 
turnaround radius and the turnaround overdensity (see points (1) 
and (2)). We say that a cluster (either regular or irregular) fits our 
model when its mass profile satisfies everywhere, in the considered 
region, the following condition (i is the cluster label, and j is the 
label of the cluster shell considered): 



logi 



.MNE(r,,,) 



0.5 ^ f, , < 2. 



(35) 



where e defines the amplitude of the agreement band, and the 
0.5 ^ fij ^ 2 interval is defined in order to reproduce the non- 
equilibrium region; its amplitude is necessarily related to the ampli- 
tude of the non-equilibrium regions of our clusters. We first adopted 

and 



the normalization in equation 1 28 s 
into equation l |35| ). We performed 1 
e = 0.10, e = 0.15, and e = 0.20; 



and introduced f, , „„_ , 
' three test for each cluster, with 
these values correspond to a 
maximum ratio between data and model of about 1.25, 1.4, and 1.6, 
respectively. The results of this comparison are reported in Table|2] 
For the three choices of e (listed in columns), the first two rows 
indicate the number of regular and irregular clusters which agree 
with our model, respectively, while the third row indicates the total 
number of clusters which agree with the model. Most clusters seem 
to be in good agreement with our model and therefore evidence the 
existence of a common profile in the external region. As expected, 
equation |9| and equation l |10^ are particularly suitable to describe 
regular profiles, while poorly fit parts of the irregular profiles. In 
fact, in most of the cases of poor fit the discordance is present only 
in the extreme outskirts of the clusters, and does not bias the esti- 
mation of the overdensity and the mass at the turnaround radius. 
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Table 3. Agreement with the model for different types of mass profiles. We 
used the condition in equation |35} with the normalization in equation \29\ . 

Agreement in the whole non-equilihrium region 





e = 0.10 


e = 0.15 


e = 0.20 


Regular 
IiTegular 
All 


23/72 (32%) 
12/42 (29%) 
35/114 (31%) 


48/72 (67%) 
21/42 (50%) 
69/114 (61%) 


64/72 (89%) 
28/42 (67%) 
92/114(81%) 


Agreement at turnaround radius rt 

6 = 0.10 e = 0.15 e = 0.20 


Regular 
Irregular 
All 


43/72 (60%) 
16/42 (38%) 
59/114 (52%) 


67/72 (93%) 
25/42 (60%) 
92/114(81%) 


71/72 (99%) 
32/42 (76%) 
103/114 (90%) 



The normalization procedure used so far implies the knowl- 
edge of the infall velocity profile of clusters in order to compute 
rt;i, Sf.i, and Mt-i- But the infall velocity profile cannot be com- 
puted directly from observations, since we know only the line-of- 
sight velocity coinponent. So, to better estimate the model reliabil- 
ity when applied to observed clusters, we should adopt the normal- 



ization in equation \29\. We therefore introduced and M^'^j 



into equation l[35j. The results of the comparison between the clus- 
ter profiles and those of our model are reported in Table[3] The over- 
all agreement between the profiles is worse than that obtained via 
the previous normalization procedure; however, despite the vari- 
ance among the infall velocity profiles of the clusters, about 80% 
of all our profiles succeed to be well described by our model in 
all the non-equilibrium region (with a maximum uncertainty corre- 
sponding to e = 0.20). 

The second panel of Table |3] displays the agreement between 
data and model at the turnaround radius rj_. In this case, we re- 
stricted the agreement interval in equation liSt to r'^' — rt.-RM- As 
one can see, the difference between regular and irregular profiles is 
smaller than it is in the previous case. Despite the uncertainty due 
to the variance among clusters, our model is able to correctly es- 
timate the cluster turnaround masses in more than 80% of cases, 
within an agreement amplitude e = 0.15. 



4 CONCLUSIONS 

We analized a large sample of simulated galaxy clusters in order to 
reconstruct the mass profile in the non-equilibrium region, where 
the galaxy dynamics is dominated by an overall infall motion to- 
wards the cluster centre. Within the assumptions of the spherical 
infall model, the turnaround overdensity St can be theoretically 
computed as a function of only the matter density parameter Qo, 
assuming a spatially flat universe. We obtained the overdensity 
St — 6 ~ 15, depending on the infall velocity profile we adopted. 

We interpolated the infall velocity profile of member galax- 
ies extracted from the simulated clusters of our catalogue, and we 
showed that: 

(i) The turnaround radius rt can be quite well approximated by 
a multiple of the virialization radius r.^: rt — 3.5r„; 

(ii) The turnaround overdensity St is consistent with the predic- 
tion of the spherical infall model, as long as the infall velocity pro- 
file is described by the Meiksin approximation ( jVillumsen & Davis] 
[T986l >. 



turnaround mass Mt and the virialization mass M^: Mt — 1.7M„. 
Moreover, Mt turns out to depend on the 3-d DM velocity disper- 
sion within the virialization core CTi,,dm approximately in the form 
of a cubic relation. 

The turnaround values can be assumed as a suitable normal- 
ization scale for the mass profiles in the non-equilibrium region of 
clusters. We showed that the normalized mass profiles are gener- 
ally consistent with a cosmic profile, which can be described (for 
0.5 <r/rt< 2) by: 



M(r) ~ Mt exp 



0.6 



1/4 



(36) 



Points (i) and (ii) are in agreement with |Vedel & Hartwick|p998) 
and |Regos & Geller| ( |1989[ > and imply a proportionality between the 



While in the inner, relaxed or almost-relaxed regions the mass 
can be considered independent on cosmological parameters, in the 
outer regions a dependence on (even if small) has to be taken, 
at least in principle, into account. 

We used a synthetic cluster, obtained by summing all cata- 
logue clusters, to determine a robust estimate of rt and St- If we 
assume this values, our model is able to predict the mass profile in 
the non-equilibrium region for about 80% of clusters. So, it is pos- 
sible to speak about a mass profile even in the region where mass 
accretion takes places along isolated radial filaments rather than in 
a spherically symmetric way. 

Our model may be useful in observational analysis in order to 
estimate the total mass of clusters using the redshift-space distribu- 
tion of galaxies. The method is the following: 

(i) one estimates the virialization radius and the virialization 
mass from the galaxy velocity dispersion in the cluster core; 

(ii) using equation l[23j and equation p.5\ one computes the 
turnaround radius rt and the turnaround mass of the cluster Mt; 

(iii) once the turnaround radius and the turnaround mass are 
known, one can estimate the mass profile in the non-equilibrium 
region using the exponential law in equation l |36| l. 

The advantage of this approach lies in the possibility to estimate 
the mass profiles up to the far outskirts of clusters, where the caus- 
tic pattern is not generally recog nizable (Diaferio & Geller 1997J 
|Diaferiop999^ ; up to now, these cluster outer volumes have been 
usually neglected in the evaluation of cluster total masses. 

Actually, the turnaround mass is a more exhaustive evaluation 
of the total mass of the cluster. The steps leading to it consist in 
the abovementioned points (i), (ii), and (iii), and are expected to be 
applied to observed clusters. 
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APPENDIX A: THE SPHERICAL INFALL MODEL IN A 
FLAT UNIVERSE 

In this appendix we briefly recall some results of the spherical infall 
model which are useful for our discussion. We consider a spherical 
density perturbation in a flat universe with cosmological constant 
(Qo + Qa = 1) and describe it as a Rriedmann universe on its 
own. Both the radius R of the perturbation and the universal scale 
factor a are normalized with ao, the present-day scale factor, and 



treated as adimensional quantities. The redshift is therefore defined 
as z — a^^ — 1. 

Our first aim is to compute the radius and the overdensity of 
the perturbation at the present day as a function of the correspond- 
ing primordial values. We use for this purpose the Rriedmann equa- 
tion for the perturbation: 

.R. (Al) 



4nG „ Ac^ 

P-RH 

3 3 



Here G is the gravitational constant, A is the cosmological con- 
stant, and p = pbg{l + S) is the density of the perturbation. Equa- 
tion \Al\ has no analytical solution, but can be solved numerically 
between an initial time ti and the present time tg, assuming the 
well-known Friedmann solution for a: 



a{t) 



no 



1 -fio 



sinh 



1/3 



(A2) 



We will hereinafter adopt the subscript in and to denote the ini- 
tial and the present value of quantities, respectively. We choose the 
initial time so as to obtain ~ ai„. In the matter-dominated era, 
we have po^io — pinO-L and af„ = pbg,o/ phg,,n- Equation i Al 
can therefore be rewritten as follows: 



F 



-2^F + 
a 



S^o - 1 F + 



r2o(i + Si, 



2a3 



-F"^ = 0, (A3) 



where F = R/a, and the dots denote first- and second-order 
derivatives with respect to r = Hot. We favour equation ^A3\ 
because it gives stabler results when integrated by computational 
means. We have Fi„ ~ 1 and Fo = Ro, which yields, in a comov- 
ing framework: 



ro = For.n, 

l + So = Fo'^il + S.n), 



(A4) 
(A5) 



where r is the comoving radius of the perturbation. 

Our second aim is to determine an analytical expression for 
the turnaround radius and the turnaround overdensity of the pertur- 
bation at the present time. In this case, we study the variation of 
R with respect to a, which can be expressed as follows l (Peebles| 
[I984| >: 



dR 
da 



R~^ +uJoR^ - 
+ u)oa? 



(A6) 



Here loo = — 1, and k > parametrizes the overdensity of 
the perturbation. When k is large enough, the perturbation expands 
until it reaches a maximum radius Rmax and then recoUapses. The 
value of Rmax comes as a solution of the following third-degree 
algebraic equation ( |Eke et al.|1996[ ): 



ujoR„ 



I^Rmax ~t~ 1 — 0, 



(A7) 



where the condition k ^ (9ti;o/4)^'''' is required for positive real 
solutions to exist. If so, we obtain: 



4k \ 



1/2 



Rmax(t^) = ( COs[6'(k)], 



where 



9{n) = - 27r — arccos 



1/2 



(A8) 



(A9) 



with 7r/2 5C ^ 27r/3. The redshift of maximum amplitude 

Zmax and the redshift of virialization of the perturbation can 
both be computed by numerical integration of equation \A6): 



10 Guido Cupani, Marino Mezzetti and Fabio Mardirossian 



1/2 

sinh[(j!>„j„j;(K)] 



,1/2 



where 



sinh[<j!>i,(K)] 



3 1/2 

2 ° 



J? 



K + 1 



1/2 



(AlO) 



(All) 



(A 12) 



and <j)v{n) ~ 2<j)max{i^), since the integral term must be taken 
twice to consider both the expansion phase and the collapse phase. 

The present-day turnaround radius rt is defined as the radius 
of a perturbation which is now reaching its maximum amplitude. 
Conversely, the present-day virialization radius is defined as the 
radius of a perturbation which reached its maximum amplitude in 
the past and is now setting to equilibrium after collapse. Let nt and 
K„ be the overdensity parameters of this two pertubations, respec- 
tively. We obtain: 



— Rmax 
1 



R max 



2 - 77„/2 

where rj^ = 2ijJoRmax{K-v)^ jLahav et ah 



(A 13) 
(A14) 



1991 



(All 



(ftt) = into equation i AlO' and 

we obtain 4>m.ax{i^t) ~ arcsinh(tt;Q''^) and 



Substituting 
into equation 



(l)v{nv)/2 = arcsinh(tJo'^^)/2. Since Rmax(rt) and 7?™a3;(r„) are 
known from equation \A&\ , we can use these relations to evaluate 
Kt and Kt,, and consequently rt and r^- We obtain in particular, as 
an original result: 

1/2 o 

(A15) 



2 — UJoRmaxiK-v)"^ 



1 — iJ->oRmax{K,y) 



3 1/2 , 
Kt, COS f 



where 9t = 6 (hit) and 6'„ 
obtain 



Since 1 + St — r^.^ , we also 



l + <5t 



4Kt (cos6't) 



3/2 



(A 16) 



APPENDIX B: INVERSION OF THE YAHIL'S FORMULA 

Let /y be the Yahil approximation function, defined in equation 
l |18^ . We consider its first-degree Taylor series expansion around 
the point 5q : 

fY{S) = h{5) + 0{5 - 5o)\ (Bl) 



A(5) = /y(5o)+ 



{5^5o). (B2) 

If the second- and higher-order terms are negligible, we can sub- 
stitute equation ( |Bl[ l into equation l[3j, thus obtaining a linear 
equation which provides 5 as a function of So and Q,Q''^''^Vr/ Hor. 
Choosing 5q = 10, the non-linear terms turn out to be negligible in 
the range where turnaround occurs: 

|C((5- 10)^1 sC 0.05/y(5), 7^5^20. (B3) 
In this case we can write: 

Sy = d{fY) ~ — 11 ' Qo — — , (B4) 

Li nor 11 

giving 5Y.t ^ 11 (when Vr = Hor) and fl'^''^ /y {Sy ,t) ^ 0.998, 
very close to unity. 



